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Abstract 

A  mixed  analytical/numerical  method  is  developed  here  to  solve  the  low  Reynolds  number  k- 
epsilon  turbulence  model.  In  this  method  the  advection-diffusion  part  is  solved  numerically, 
while  the  source  terms  are  split  into  two  parts:  one  part  is  solved  analytically  and  the  next  is 
solved  numerically. 
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1.  Introduction 

Because  of  the  recent  enormously  progress  in  the  capability  of  computers,  low  Reynolds  two- 
equation  turbulent  models  become  more  and  more  welcome  in  engineering  fluid  computation. 
In  the  past,  various  forms  of  low-Reynolds-number  k  -  s  turbulent  models  have  been  proposed. 
The  detail  of  two-equation  models  and  low  Reynolds  corrections  are  presented  in  Chapter  14 
(Turbulence  Modeling  and  Simulation)  of  handbook  [1], 

Though  mathematically  the  k  —  e  model  is  well-posed[2],  the  strong  nonlinearities  may 
interact  with  numerical  errors  in  such  a  way  that  computation  may  break  down  easily..  A 
typical  behavior  of  unstable  computations  involves  the  loss  of  positivity  of  k  or  e,  though  the 
original  differential  equations  have  positive  solution  [3].  The  appearance  of  negative  values 
changes  the  sign  of  several  terms  in  the  models,  so  that  turbulent  quantities  may  increase 
unboundedly  [4],  Even  though  the  two  turbulence  equations  can  be  solved  exactly  in  the  same 
manner  as  the  mean  flow  equations,  it  has  been  found  that  such  a  method  often  leads  to  an 
unstable  solution,  or  even  incorrect  solutions  [5].  The  damping  functions  in  the  low  Reynolds 
turbulent  models  improve  the  model  prediction  capability  for  near  wall  flow,  but  also 
introduce  more  severe  numerical  stiffness  for  the  source  terms. 

There  are  a  huge  amount  of  numerical  methods  for  compressible  Navier-Stokes  equations 
coupled  with  two  equation  models  [6-9].  The  two-equation  turbulent  model  is  a  typical 
example  of  partial  differential  equations  with  source  terms.  Great  progress  has  been  made  in 
efficient  treatment  of  the  source  terms  [10-12].  Helzel,  LeVeque,  and  Wamecke  [12]  treated 
chemical  reacting  flow  with  an  Arrhenius  law  for  the  source  term  by  mixed  method  in 
detonation  waves  computation.  In  [13],  a  mixed  analytical/numerical  method  for  oscillating 
source  terms  has  studied.  In  this  method  the  advection-diffusion  part  and  the  source  terms  are 
treated  separately  through  operator  splitting.  The  advection-diffusion  part  (PDE)  is  integrated 
numerically  while  the  source  term  part  (ODE)  is  integrated  analytically.  Hence  this  method  is 
called  mixed  analytical/numerical  method.  The  mixed  method  performs  well  for  partial 
differential  equations  with  source  terms,  in  which  the  time  scale  of  source  term  (S-scale, 
denoted  Ts)  is  much  smaller  than  the  mean  flow  scale  (M-Scale,  denoted  Tty)  inherent,  to  the 
advection-diffusion  part.  Furthermore,  the  mixed  method  has  been  extended  to  the  implicit 
solution  of  high  Reynolds  number  and  compressible  turbulent  flows  [14].  Numerical  results 
show  the  mixed  method  can  give  robust,  steady  and  fast  convergent  solution. 

In  this  paper,  we  extend  the  mixed  method  to  low  Reynold  number  turbulent  models.  The 
essential  new  feature  of  the  mixed  method  for  low  Reynolds  number  turbulent  models  lies  in 
the  treatment  of  the  source  terms  which  contain  new  damping  functions.  With  respect  to  the 
high  Reynolds  number  counterpart,  the  low  Reynolds  number  turbulent  models  retain  the 
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orginal  source  terms  modified  with  damping  factors  and  sometimes  contain  additional 
damping  terms.  In  the  mixed  method  proposed  in  this  paper,  the  damping  factors  are  treated  as 
constant  at  each  time  step,  so  that  the  main  part  of  the  source  terms  are  still  analytically 
integrable  at  each  time  step.  The  additional  damping  terms  are  treated  numerically. 

2.  Governing  Equations 

The  governing  equations  are  obtained  by  Favre  Averaging  the  Navier-Stokes  equations 
and  modeling  the  Reynolds  stress.  In  conservative  form  these  equations  are  written  as 


dU  dFc  dGc 

- +  — -  + — - 

dt  dx  dy 


dF ,  8G„ 

-  +  - 


dx  dy 

In  this  paper,  we  use  the  Hwang-Lin  low-Reynolds  number  k-  &  turbulence  model1 14  as 
example. 
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The  production  term  Pk  is  given  by 
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The  eddy  viscosity  is  calculated  as 
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The  functions  II  and  are  defined  as 
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The  coefficients  and  damping  functions  are 

Cm=0.09,  Ce1=1.44,  Ce2=1.92 
4=l-exp[-0.01v,-0.008vi] 

crk  =  1  -  exp[-0.0 1  -  0.008j]  ] ,  <re  =1.3-1.0exp[-^/10] 

y; =i.o,  f2= i.o 


&2v(afj 


£  = 

where  yx  = 


vk/S'c 

3.  Numerical  Method 

The  two  set  equations  (1)  and  (2)  are  solved  separately.  The  convective  numerical  flux  at  the 
cell  interface  is  evaluated  using  Roe’s  approximated  Riemann  solver  with  MUSCL  treatment 
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to  achieve  second  order  accuracy.  The  simplified  multistage  scheme  with  four  stages  is 
applied  to  advance  in  time. 

The  advection-diffusion  part  of  the  two-equation  model  is  discretized  in  the  same  way  as 
for  the  Navier-Stokes  equations.  Only  the  source  term  needs  special  treatment.  The  numerical 
treatment  of  source  terms  of  the  turbulence  model  equations  is  of  mayor  importance  for  the 
stability  of  this  scheme  [16].  The  fundamental  means  is  to  treat  the  negative  source  terms 
implicitly  and  the  positive  terms  explicitly[7].  However,  the  time  step  still  has  to  be  small 
enough  in  order  to  obtain  realistic  values  of  k,  s  and  pt.  According  to  the  different  treatment  for 
Jacobian  matrix  for  implicit  part,  there  are  three  method:  point  implicit  method[6], 
approximate  Jacobian  method[5,17]  and  exact  Jacobian  method[18].  In  this  paper,  we  adopt 
the  second  method. 

4.  Mixed  Analytical/numerical  Method  for  low  Reynolds  number  turbulent  models 

The  turbulence  model  equations  contain  advection-diffusion  operators  and  source  terms.  In  the 
mixed  method  the  advection  diffusion  part  (PDE)  is  integrated  numerically,  while  the  source 
term  (ODE)  is  integrated  analytically. 

4.1  Treatment  of  the  source  terms 

For  the  case  of  low  Reynolds  number  models  it  is  not  always  possible  to  integrate  the  source 
terms  analytically.  With  respect  to  the  high  Reynolds  number  counterpart,  the  low  Reynolds 
number  turbulent  models  retain  the  original  source  terms  modified  with  damping  factors  and 
sometimes  contain  additional  damping  terms.  In  the  mixed  method  proposed  in  this  paper,  the 
damping  factors  are  treated  as  constant  at  each  time  step,  so  that  the  main  part  of  the  source 
terms  are  still  analytically  integrable  at  each  time  step.  The  additional  damping  terms  are 
treated  numerically. 

Following  the  original  construction  of  the  mixed  method,  we  first  consider  the  ODE  due 
to  the  source  terms 


dUT 

dt 


=  Ss(U)+SD(U) 


(V) 


Here  Ss(U)  is  the  standard  part  which  is  similar  to  the  source  term  of  the  high  Reynolds 
number  model,  modified  only  by  adding  damping  factors,  and  SD(U)  stands  for  the  additional 
damping  terms.  Precisely,  Ss( U)  is  defined  by 


Ss(U)  = 


Pk~> 


(8) 


yfCel  k  Pk  f2Ce2  k  j 

The  additional  damping  terms  SD(U)  are  treated  numerically,  while  the  standard  part  Ss(U) 
is  solved  analytically  by  considering  the  following  ODE 
d  pk 


dt 

dpffc 
dt 

with  the  parameters 


=  pt-pffc 


(9) 
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treated  as  constant  at  each  time  step.  Hence  the  standard  part  should  have  the  same  form  of 
analytical  solution  as  the  high  Reynolds  number  model[13],  and  this  analytical  solution  for  (16) 
can  be  written  as 


pk  =  K((pk\,(p^,t)  =  (pk\2aFr 
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Here  ko  and  Sfc-  stand  for  initial  value  of  turbulent  kinetic  energy  k  and  dissipation  rate  s. 

4.2  Mixed  method 

Now  the  mixed  method  for  low  Reynolds  number  models  can  be  written  as 

[/-A#.[^]"]A(^)’=At./j”+Ss((^),,)  +  SD((^)”)  (12) 

where  s  ^;jl  y  j  is  due  to  the  additional  damping  terms,  and  <r  is  the  analytical  solution 

of  the  ODE  due  to  the  standard  source  terms.  This  analytical  solution,  when  limited,  can  be 


written  as 


u1 


minj K((pk)"j,(p£ty'j,At')-(pk)tj,  8 
min 


(13) 


The  parameter  8  limits  the  maximum  of  the  source  term  integral.  The  choice  of  8=0.1  works 
well  in  our  numerical  experiments. 

5.  Numerical  Results 

5.1  Flat  plate 

We  first  calculate  a  turbulent  flow  past  a  flat  plane.  The  grid  is  80x41.  The  distance  of  the  first 
grid  next  to  the  wall  is  no  more  than  y+=0.45.  The  Reynolds  number  is  Re=1.0xl07,  the  inflow 
Mach  number  is  Ma=0.2.  The  convergence  history  of  mixed  method  is  compared  with  that  of 
traditional  method  in  Fig.l.  The  latter  can't  converge  with  the  same  condition.  The  inner  layer 
velocity  profile  of  flat  plane  flow  is  shown  in  Fig.2,  comparing  the  van  Driest' s  and  Spalding's 
theory.  The  figure  gives  the  velocity  profile  of  7  sites  in  x  direction.  The  agreement  is  also 
good. 

5.2  Transonic  Diffuser 

Consider  transonic  flow  with  a  weak  shock  through  a  converging  diverging  diffused  19].  This 
configuration  has  an  entrance  to  throat  area  ratio  of  1.4,  an  exit  to  throat  area  ratio  of  1.5.  The 
corresponding  Reynolds  number  is  9.370xl05  and  the  Mach  number  is  Min=0.9.  The  total 
pressure  at  inflow  is  1.349xl05Pa,  the  static  pressure  at  the  outflow  is  l.llxl05Pa.  These 
flows  were  characterized  by  the  ratio,  R,  of  exit  static  to  inflow  total  pressure.  For  the  weak- 
shock  case  the  value  of  R  was  0.82. 
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The  convergence  history  is  shown  in  Fig.3,  comparing  that  with  traditional  method.  We  see 
the  computation  with  traditional  method  fails  to  converge,  while  the  mixed  method  not  only 
converge  steadily,  but  also  converge  faster  than  the  traditional  method. 


Fig-1  Fig.2 

Fig.  1.  Comparative  convergence  history  for  flat  plane 
Fig.2.  Velocity  profile  of  flat  plate  flow 


Traditional 
Mixed  Meth 


Vi  | 


0  5000  10000  15000 

Iterative  Time 

Fig.3 

Fig.3.  Comparative  convergence  history  for  transonic  diffuser 
Fig.4  Velocity  profiles  at  four  axial  locations  for  transonic  diffuser 

5.3  RAE2822  Airfoil 

The  RAE2822  airfoil  has  been  extensively  used  for  the  validation  of  Navier-Stokes  codes 
applied  to  transonic  flow.  The  key  flow  condition  for  this  test  case  is:  Free  stream  Mach 
Number  Maoo=0.73,  Reynolds  number  Re=6.5xl06,  and  angle  of  attack  in  degree  a=3.19.  The 
convergence  histories  with  mixed  method  and  with  conventional  numerical  method  are 
compared  in  Fig.5.  The  first  one  converges  faster.  The  pressure  coefficient  and  skin  friction 
coefficient  with  experimental  results  are  shown  in  Fig.6  and  Fig.7  respectively,  comparing  the 
conventional  method  results.  The  numerical  results  with  mixed  method  agree  with  experiment 
data. 


Fig.  5.  Comparative  convergence  history  for  RAE2822  airfoil 
Fig.  6.  The  pressure  coefficient  profiles  of  RAE2822  airfoil 
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Fig.  7.  The  skin  friction  coefficient  profiles  of  RAE2822  airfoil 
5.4  NLR7301  Wing-flap  Airfoil 

The  flow  around  the  two-element  NLR7301  airflow  with  2.6%  gap  has  been  computed  for  the 
flow  condition:  Re=2. 51x106  and  angle  of  attack  a=13.1.  Fig.8  gives  the  pressure  distribution 
on  wing  and  flap. 


Fig.  8.  The  pressure  coefficient  profdes  of  NLR7301  airfoil 

6.  Concluding  Remarks 

The  mixed  analytical/numerical  method  has  been  extended  to  the  numerical  solutions  of  low 

Reynolds  number  k-s  turbulence  models.  The  mixed  method  is  applied  to  FIwang-Lin  low 

Reynolds  turbulent  model  and  several  test  problems.  The  numerical  results  show  that  the 

mixed  method  is  numerically  more  robust  than  the  traditional  pure  numerical  method. 
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